N90- 19407 


1989 

NASA/ASEE SUMMER FACULTY FELLOWSHIP PROGRAM 

MARSHALL SPACE FLIGHT CENTER 
THE UNIVERSITY OF ALABAMA IN HUNTSVILLE 

EVALUATION OF THE HEAT TRANSFER MODULE (FAHT) OF FAILURE 
ANALYSIS NONLINEAR THERMAL AND STRUCTURAL INTEGRATED 

CODE (FANTASTIC) 


Prepared by: 

Majid Keyhani 

Academic Rank: 

Assistant Professor 

University and : 
Department 

University of Tennessee, Knoxville 
Mechanical & Aerospace Engineering 

NASA/MSFC: 


Laboratory: 

Division: 

Branch: 

Structures and Dynamics 
Thermal Engineering & Life Support 
Thermal Analysis 

MSEC Colleague: 

Kenneth E. McCoy 

Date: 

September 12, 1989 

Contract No. : 

The University of Alabama in 
Huntsville, 

NGT-01-008-021 


xvn 




Evaluation of the Heat Transfer Module (FAHT) of Failure Analysis Nonlinear 

Thermal and Structural Integrated 
Code (FANTASTIC) 

Majid Keyhani 

Assistant Professor of Mechanical and Aerospace Engineering 
University of Tennessee 
Knoxville, Tennessee 

ABSTRACT 

The heat transfer module of FANTASTIC Code (FAHT) is studied and evaluated to 
the extend possible during the ten weeks duration of this project. A brief background 
of the previous studies is given and the governing equations as modeled in FAHT are 
discussed. FAHT’s capabilities and limitations based on these equations and its 
coding methodology are explained in detail. It is established that with improper 
choice of element size and time step FAHT’s temperature field prediction at some 
nodes will be below the initial condition. The source of this unrealistic temperature 
prediction is identified and a procedure is proposed for avoiding this phenomenon. It 
is further shown that the proposed procedure will converge to an accurate prediction 
upon mesh refinement. Unfortunately due to lack of time, FAHT’s ability to 
accurately account for pyrolysis and surface ablation has not been verified. 
Therefore, at the present time it can be stated with confidence that FAHT can 
accurately predict the temperature field for a transient multi-dimensional, 
orthotropic material with directional dependence, variable property, with nonlinear 
boundary condition. Such a prediction will provide an upper limit for the 
temperature field in an ablating decomposing nozzle liner. The pore pressure field, 
however, will not be known. 
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Permeability matrix 

Components of the Permeability matrix 

Defined by Eq. (FAHT-5) 

Specific heat 

Thermal Conductivity 

Thermal Conductivity matrix 

Components of the Thermal Conductivity matrix 

Pore pressure 

Volumetric heat source/sink 

Heat convected by pyrolysis gas, defined by Eq. (FAHT-2) 
Specific gas constant 
Rate of surface recession 
Time 

Time step size in an analysis 
Temperature 
Gas Temperature 
Initial Temperature 
Surface Temperature 
Velocity vector 

Vector differential operator del 

Coordinate in 3-D space 

A generic cartesian coordinate 

Element thickness in the direction of heat flow 

Thermal diffusivity 

Penetration depth 

Defined by Eq. (FAHT-3) 

Fmissivity 

Eigen value of the exact solution 
Dynamic viscosity 
Pi 

Density 

Porosity 

Volumetric mass source/sink 
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I. INTRODUCTION 


1.1 Overview 


The FANTASTIC code is in its final stage of development by Failure Analysis 
Associates (FaAA) for Marshall Space Flight Center. The code is presently 
under review and evaluation for its capabilities. The intended use of the Code 
is to increase the capabilities and accuracy of the thermal and structural 
analysis of solid rocket motor nozzles. The code consists of three modules for 
thermochemical analysis, heat transfer and mass diffusion analysis, and 
structural analysis. 

11.2 Objectiv es of the Present Work 

The present work is limited to the evaluation of the heat transfer and mass 
diffusion module (FAHT) of the FANTASTIC Code. In order to proceed with 
the stated task within the rather short time period of the project, the following 

sub tasks were chosen: 

1 . A brief review of previous efforts in the area of code development. 

2. Verification of the various required capabilities of the FAHT module, 
such as transient, nonlinear boundary condition and variable property 
solution routines. 

3. A limited attempt at use of FAHT for prediction of the temperature 
distribution in MNASA nozzle. 
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H. BACKGROUND 


Charring ablators have proven to be one of the most successful thermal protection 
systems for applications with high thermal loading such as reentry and solid rocket 
nozzles. The materials used are a combination of plastics that decompose to a porous 
char zone and pyrolysis gases. The pyrolysis gases transport energy to the surface 
under thermal loading thereby reduce the rate of energy input into the virgin 
material. The decomposition process along with endothermic reactions of the 
pyrolysis gases with the carbon in the char zone are of further help in reducing the 
rate of energy input into the virgin material. A cross-sectional view of a charring 
and decomposing ablator of a nozzle is depicted in Fig. 1. The events taking place in 
the char and pyrolysis zones can be summarized as follows: 

Char Zone 


a. Thickness of this zone is about 2 to 4 mm depending on the material used. 

b. 1 yrolysis gases flowing through this porous zone have a cooling effect due to 
convective transport of energy to the surface. 

c. Pyrolysis gases are not frozen in this region and endothermic reactions with 
carbon of the char zone takes place, resulting in a "cooling” of the char zone. 

d. Chemical reactions between the char zone and the pyrolysis gases result in a 
constant change in the porosity and the permeability of the char layer. 

Pyrolysis Zone 


a. I hickness of this zone in typical ablators is about 1 to 2 mm. 

b. Porosity and permeability in this region are changing rapidly. 

c. A volumetric energy source/sink is present due to decomposition of the virgin 
material. 
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Chemical species 

adiation diffusion 
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Figure 1. Cross-sectional view of a charring and decomposing ablator of a nozzle. 



The above summary of the phenomena taking place in a charring ablator clearly 
shows the complexity of the problem. Other factors that add to the difficulty of 
accurate prediction of temperature and pore pressure values in a charring ablator 
are: 

a. Temperature variation in space and time is substantial, therefore the 
evaluation of the thermophysical properties at the proper temperature is very 
important. 

b. The material is orthotropic, therefore the conservation equations for energy 
and momentum are more difficult to model. 

c. As mentioned earlier porosity and permeability are changing with time. 
Therefore, they are variables in the momentum equation. The porosity 
variation can be accounted for via its relation to density variation governed by 
the rate equation of the Arrhenius form. The modeling of the variable 
permeability for an orthotropic material, however, is not trivial. It appears 
that no empirical relation for relating the permeability of this class of 
materials to its porosity is available at the present time. 

d. The nonequilibrium chemical reaction between the pyrolysis gases and the 
char zone has not been throughly investigated for all the candidate materials 
for nozzle liner. 
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11.1 Previous Work 


A review of the literature indicates that research and code development effort in the 
area of thermal analysis of charring and decomposing ablators had a rather high 
priority for NASA during the decade of 1962-1972. It appears that the development 
of new candidate materials for a charring-decomposing ablator may require 
initiation of a rigorous research program in this important area. A brief summary of 
the previous work is given in this report. 

• Acrotherm charring material thermal response and ablation program (CMA, 
Aerotherm report 75 148, 1975) 

One-dimensional transient finite difference model 

No resistance to flow of pyrolysis gases 

Frozen flow of pyrolysis gases through the char zone 

• NASA report NASA TN-D-6895 (1972) 

Two-dimensional transient axisymmetric finite difference model 
Resistance to flow of pyrolysis gases incorporated via Darcy’s law 
Frozen flow of pyrolysis gases through char zone 

• NASA reports NASA TN D-1370 (1962) CHAP I and NASA TN D-2976(1965) 

ciiAPir 

One-dimensional Transient finite difference program 
No resistance to flow of pyrolysis gases 

Frozen Flow of Pyrolysis gases through the char zone (CHAP I) 

- Chemical reaction between the pyrolysis gases and the char incorporated 
(CHAP II) 

• NASA report NASA TN I) 6085 (1970), COSMIC 

One dimensional transient finite difference model 

No resistance to flow of pyrolysis gases 

Frozen flow of pyrolysis gases through the char zone 
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• NASA report NASA CIM 903 (1971) 

- One dimensional steady-state finite difference model 

- Modified form of Darcy’s equation used to model resistance 
to flow of pyrolysis gases 

1 low of pyrolysis gases through the char zone analyzed as 

1. Frozen 

2. Equilibrium 

3. Non-Equilibrium 

H.2 Governing Equations 

In this section the energy, momentum and mass conservation equations as they 
appear in the FANTASTIC/FAHT theoretical manual (version 1.0) will be reported 
and discussed. 

Energy equation: 


PC 


p 


<rr 

fit 



(FAHT-1) 


The volumetric source/sink term (Q) in Eq. (FAHT-1) accounts for the energy 
associate with the decomposition process. It should be noted that the energy 
transport term due to convection of pyrolysis gases does not appear in Eq. (FAHT-1). 
P All T theoi ectical manual states that the heat convected by pyrolysis gases is 
calculated by 


Q c = P C P v ‘ VT, (FAHT-2) 

and can be accounted for via: NONLINEARHEATBC, THE RMALG APCONT ACT 
OR LU MPE DHE ATC APACJTY Options. 
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Proper form of the energy equation is given by 


nr a / aT \ 

" IT +pC -* v ' VT= — ( K , — ) + Q 0) 

L dx ] \ aXj / 

It is very important to note that the temperature pressure and velocity in a 
decomposing ablator are undergoing substantial changes with a rather small change 
in time. Literature contains a very large volume of papers dealing with various 
convection heat transfer problems. This author is not aware of a single case where 
the convection term is not explicitly treated in the energy equation as is the case 
with Eq. (FAHT-1). 


In verbal communication with Dr. McCoy (NASA/MSFC-ED64) FaAA Personnel 
have stated that the heat convected by pyrolysis gases (Eq. FAHT-2) is accounted for 
as part of the load term Q. In that case the energy equation is not related accurately. 


Momentum Equation: 

Y 

V = - B V P 1 

4>p 


(FAHT-3) 


dP 

Y - K , when ideal gas or y = 

dp 

lhe momentum equation as given by Eq. (FAHT-3) is dimensionally inconsistent. 
The proper form of the momentum equation based on the Darcy law is 

1 

V = B • (VP) (2) 

4> p 

Mass Conservation Equation: 
d 

[4)p] +V-($pV) = $ (FAHT-4) 
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In order to obtain an explicit equation for pore pressure, FAHT theoritcal manual 
states that Eqs. (FA1IT-3) and (FAHT-4) are combined to get. 


d 

— [C m PJ = 
at 

a 

/ B tJ dP 

dXj 

' y dXj 

where C„, = $ 

dp 

4> 

nr — 

dP 

RT 


) 


+ <t> 


(FAHT-5) 


It should be noted that the pore pressure equation as given by Eq. (FAHT-5) is also 
dimensionally inconsistent. The proper form of the pore pressure equation is 


dt dx 

I 



(3) 


P 

where p = 

R T 

g 


The problem with Eqs. (FAHT-3) and (FAHT-5) may be due to a typographical error, 
and they are correct in the program. If that is the case, then it is an indication of a 
poor editing job of various versions of manual up to version 1.0 which have been 
released by FaAA. 


TT.3 FAHT’S Capabilities and Limitations 

In order to get an estimate of FAHT’s potential capabilities it is necessary to come to 
a concl usion about the governing equations as described in FAHT’s theoretical 
manual (version 1.0). The following assumptions are made with regards to modeling 
of the governing equations in the FAHT module. 

(a) 1 he convected energy by the pyrolysis gases are accounted for as part of 
the load term (Q) in Eq. (FAHT-1) 

(b) The problem with Eqs. (FAHT-3) and (FAHT-5) is typographical, and these 
equations are indeed dimensionally consistent in the program. 

Eased on the above assumptions, FAHT’s potential capabilities in analysis of a 
charring decomposing ablator are as follows: 
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(i) The temperature field solution will have some errors due to lack of explicit 
treatment of the convection term. 

(ii) The pore pressure equation, as given by Eq. (3) is nonlinear. This 
equation, however, is treated as a linear equation in FAHT. Therefore, the 
pore pressure solution will not be accurate, which will result in an 
inaccurate velocity solution. The inaccuracy associated with the velocity 
solution, in turn, will increase the error in the temperature solution. 

(iii) The flow of pyrolysis gases through the char zone are modeled as a non- 
reacting low (frozen). Therefore, the endothermic reactions which take 
place in this zone can not be accounted for. Furthermore, the changes in 
the porosity and permeability in this zone cannot be calculated and 
accounted for. 

(iv) The permeability of the material is assumed to have the same value in the 
virgin, decomposition and char zones. As mentioned earlier, the 
permeability is changing rapidly with time in the decomposition zone. 
Nonetheless, FAHT can not account for this variation. 

(v) There are no provisions for accounting for the initial porosity of the virgin 
material. 

(vi) The momentum equation is based on the Darcy law. However, it is known 
(NASA CR- 1 903, 197 1) that the inertial effects play an important role due 
to relatively high mass fluxes of the degradation products. An accurate 
modeling of the momentum equation requires the inclusion of the 
Forchheimer term. Another motivation for the use of Forchheimer- 
extended Darcy equation of motion for flow through porous media is the 
following. One of the objectives of the exploratory test program of the Solid 
Propulsion Integrity Program (SPIP) is to provide empirical relations for 
permeability of various candidate materials. There may be cases where 
the permeability should be determined via the Forchheimer-extended 
Darcy equation. In that case we do not have any analysis tool that can 
properly use the permeability data. 
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III. VERIFICATION OF FAHTS’ CAPABILITIES 


I II. 1 Transient Solution Routine 

In the process of assessment and verification of the heat transfer module of the 
FANTASTIC code, many attempts have been made at obtaining solutions for simple 
cases of transient conduction heat transfer. It appears that in many cases the 
predicted short-time temperatures at some nodes were below the initial input values. 
Thus, the code indicated that those nodes were cooling while the physics of the 
problem imposed a rather high heating rate for the material. A number of 
explanations along with corrective remedies have been forwarded by FaAA. It has 
been suggested by FaAA that for high heating rate cases (high thermal gradients) 
the input values of time step and element size should be such that the value of the 
element Fourier number is rather low. The element Fourier number is defined as 
uAL/Ax 2 . The parameter u is the thermal diffusivity of the material and can not be 
arbitrary changed. The only alternative for reducing the element Fourier number is 
then using a smaller time step and/or larger elements. The choice of a smaller time 
step and/or larger element size does not resolved the problem and indeed results in 
even lower values. Therefore, we are in a situation where we cannot get realistic 
solutions to transient problems with high heating rates. It may be noted that our 
i ntended application is precisely what FAHT apparently cannot solve. 

The source of the problem with the FAHT’s inability for solving transient problems 
with high gradients is not a "bug.” It is a rather fundamental problem related to the 
physics of the process. Before we proceed further, the terms "unrealistic” and 
"realistic” should be defined. First, note that in solving "real world” problems we do 
not know the exact answer. Thus, a realistic answer is one which makes physical 
sense to the analyst but may not be accurate. Realistic answers can be very 
dangerous. Shortly, I will explain why I consider them dangerous to an 
inexperienced analyst. 1 he unrealistic answer is very easy to detect such as the case 
of cooling nodes predicted by FAHT where the model has specified a high heating 
rate. 

Finite difference computer codes rarely give an unrealistic answer to a high heating 
rate transient problem. However, finite element codes, depending on the choice of 


XVH-10 



shape function, are prone to give unrealistic answers when the thermal gradients 
are high. As stated earlier, the problem with this class of finite element programs 
(including 1 AHT) is one of fundamental nature. A high thermal gradient results in 
a propagation of a thermal front into an isothermal domain (initial condition). For a 
given time step, this thermal front moves into the domain by a distance called 
penetration depth (6). Figure 2 depicts the movement of penetration depth into a 
domain. If 8 is smaller than the element length (Ax), then depending on the choice of 
the shape function, the nodal temperatures are determined with unrealistic values 
as depicted in Fig. 3. For further detail about this phenomena please refer to Hogge 
and Gerrekens (1982). The remedy for this problem is a choice of a larger time step 
and/or smaller elements. Obviously, if the solution routine is explicit, then the 
chosen time step should not violate the stability criterion. 

I o illustrate the fact that the suggested remedy does indeed result in a realistic 
answer by 1 AHT , let s consider a simple test problem as shown in Fig. 4. 


Test Problem 


1 ransient one-dimensional conduction heat transfer 
in a plane wall with constant properties. The 
parameters of the problem are: 

width of the wall — 0.2m 
density = 2600 kg/m 3 
specific heat = 808 J/kg-K 
thermal conductivity = 3.98W/m-K 

Boundary Conditions: 

T(x = 0,t) = T s = 3600 K 

tl'l /ilX (x = L, t), insulated back wall 

Initial Condition: 

T(x,t = 0) =T; = 300 K 

? he temperature distribution in the wall at time — 10 seconds will be presented in 
the following tables. All solutions are obtained via the implicit formulation. 
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T (0,t) = T ; 
T (x,o) = Tj 



Figure 2. Short-time transient temperature profile in a domain. 
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Input Data: 

T = 3600 K, T. = 300 K 
K = 3.98 W/m-K 
p = 2600 Kg/m 3 
C () = 808 J/Kg -K 
L = 0.2 m 
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Table 1 : FAHT temperature predictions (T (x,t) - T ( ) for various element sizes, and 


time steps at time = 10 seconds. 


Run 

No. 

X = 0.01 
(m) 

0.02 

0.03 

0.04 

No. of 
Elements 

Time Step 
(Seconds) 

FE-1 

189.4 

-120.1 

24.6 

-2.4 

20 

2 

FE-2 

261.1 

-124.3 

20.3 

-0.4 

20 

0.2 

FE-3 

233.0 

-6.9 

-0.1 

0.03 

r 40 

2 

FE-4 

272 

1.0 

0 

0 

80 

2 


I he following observations can be made with regards to the results in Table 1: 


a. A comparison of the results of runs FE-1 and FE-2 shows that when 

unrealistic (negative) nodal values are obtained, a decrease in time step does 
not help. 


b. A comparison of Ihe results of runs FE-1, FE-3, and FE-4 indicates that when 
unrealistic nodal values are obtained, increasing the number of elements 
theieby decreasing Ax yields a realistic solution. 


c. Consider that the exact answer is unknown. Furthermore, consider that the 
first choice of an analyst was the input values of run FE-4. The results of this 
run are indeed realistic and make sense. A 20cm wall is divided into 80 
elements (Ax = 2.5 mm). Moreover, a time step of 2 seconds for transition 
from zero time to 1 0 seconds is rather a reasonable choice. 


d. 


The above discussion indicates that the analyst may well choose to accept the 
results of run FE-4 as reasonable. However, the "realistic” answer of the run 
T* Iv4 is in substantial error. Indeed, the values at x - 0.01 and 0.02 are, 
respectively, 21% and 74% lower than the exact answer. This is a good 
example of the case were "realistic” answers should be treated with caution. 
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1 hope that my earlier comment about the potential danger of a "realistic” answer is 
clarified. It is reasonable to ask then, how should an analyst ensure that the 
realistic answer is reasonably accurate when the accurate answer is unknown? 

The penetration depth (8) of a thermal front into an isothermal semi-infinite domain 
can be approximated as (Ozisik, 1980) 


6 = \/l2 at (4) 

In order to obtain a realistic answer, 8 for the first time step should be large enough 
to cover a number of elements i.e., 


8 (A/) = \/ 1 2 a A / > Ax (5) 

The stability criterion for the explicit solution scheme is given by, 

Ax > \/2aA/ (6) 

Equation (6) describes an inter dependence between the time step and the element 
size. In implicit solution routines this equation is irrelevant in terms of the stability. 
However, any stable implicit solution is not necessarily an accurate solution. In 
order to ensure that an implicit solution is a fairly accurate one, gross violations of 
equation (6) should be avoided. 

The following relation which satisfies equation (5) and does not result in a gross 
violation of equation (6) is proposed 


v/aAl <Ax< v/3o A/ (7) 

The upper limit in the equation (7) ensures that 8 for the first time step covers more 
than two elements . The lower limit on Ax as given by equation (6) is relaxed by a 
factor of /2? 
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A Suggested Procedure for Obtaining Accurate Transient Results from FAHT: 


1. Start with a reasonable time step. 

2 . Select Ax by calculating the acceptable limits, and obtain the code 
prediction. 


3. 


4. 


\/ a At < Ax < V 3a At 

Reduce Ax by a factor of two and select At using the 
obtain the solution for the new Ax and At. 

Ax 2 Ax 2 

< At < 


( 8 ) 

following relation, and 

(9) 


3a a 

Compare the two results. If the changes in the temperature field are more 
than an "acceptable variation” then repeat step 3. Obviously, the degree of 
the accuracy depends on the selected criterion for the "acceptable 
variation.” 


It should be noted that the above procedure for selection of At violates the stability 
ci iterion for explicit scheme and should not be used for explicit solution. 

Lets apply the above procedure to the example problem: 

Step I : At = 2 seconds (the same as run FE-1). 

Step 2. 1.95mm < Ax < 3.37mm. Select Ax = 2.5mm (80 elements). This 

step clearly shows that for a time step of At = 2 seconds, 20 and 40 
elements runs were inappropriate. Note that the results for Ax = 
2.5 mm and At = 2 sec. are realistic (please see Table 2). 

Step 3: Double the number of elements. 

160 elements; Ax = 1.25 mm. 

0.27 < At < 0.82; choose At = 0.4 seconds. 

Step 4: A substantial change in the temperature field is observed. Repeat 

step 3. 

320 elements; Ax = 0.625 mm. 

0.07 < At <' 0.21; choose At = 0.1 seconds. 
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The results of these numerical experiments are tabulated in Table 2. The exact 
solution is also given in this table. A finite difference code based on the implicit 
formulation was programmed and the above outline for a "search” for an "accurate” 
answer was followed. The results of this program are also included in Table 2. 

TABLE 2. "Search” for Accurate Temperature Prediction, (T (x,t) - T), at Time = 


10 seconds. 


Solution 

Method 

x = 0.01 
(m) 

0.02 

0.03 

0.04 

No. of 
Elements 

Time Step 
(Seconds) 

Exact 

344.1 

3.8 

0.01 

~ 0 

- 

- 

FAHT 







FE-4 

272 

1.0 

0 

0 

80 

2 

P'E-5 

329.2 

2.75 

0 

0 

160 

0.4 

FE-6 

340.4 

3.5 

0 

0 

320 

0.1 

Finite 

Difference 







FD-1 

352.6 

14.2 

2.0 

0 

80 

2.0 

FI)- 2 

345.6 

6.1 

0.03 

0 

160 

0.4 

FD-3 

344.4 

4.4 

0 

0 

320 

0.1 

SJNDA 87 

344.9 

4.4 

0 

0 

320 

0.1 


The f A1IT results in Table 2 show that the above procedure in 3 iterations resulted 
in a fairly accurate answer (compare runs FE-6 and FD-3 with the exact solution). 
The exact solution is given in Fig. 4. 

A reflection on the number of elements or nodes used to obtain the results given in 
Table 2, poi nts to the gross inefficiency in the grid selection. For example, 320 
elements are uniformly distributed in the domain of x = 0 to 0.2. However, the 
results at time = 10 seconds show that the heat transfer is taking place in the region 
of x — 0 to 0.03, which contains 48 elements. Therefore, the remainder of the 
elements (272) are practically irrelevant. Of course, a number of techniques, such as 
deforming grid formulation, are available for efficient solution of this class of 
problems (ilogge and Gerrekens, 1982). 
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HI-2 Nonlinear Boundary Condition Routine 


In order to verify the nonlinear boundary condition routine of FAHT, the boundary 
condition of the test problem at the x = 0 surface is changed. There are various ways 
to impose a nonlinear boundary condition at this surface. The most logical one, 
considering the intended usage of the FAHT program, is a radiation boundary 
condition. I hus, the boundary condition at this surface is changed to radiation heat 
transfer from a gas at a temperature of T g = 3600 K. The surface emissivity e is 
taken to be 0.8. In section III. J it was determined that for the test problem the 
values of Ax — 0.625mm and At — 0.1 sec will result in an accurate solution by 
FAHT. Therefore, the radiation boundary condition problem is solved by FAHT with 
Hie same values for Ax and At . The predicted temperatures, T(x,t) - Ti, at time = 2 
sec for nodes from x = 0 to 2.5mm are tabulated in Table 3. In order to verify the 
accuracy of FAHT’s prediction, the predicted values obtained from SINDA-87 
program (for the same parameters) are also given in Table 3. 


The maximum variation between the results of FAHT and SINDA-87 is less than 
0.1%. It should be noted that SINDA-87 is a widely used program and its accuracy in 
solving rather complex problems has been established over the years. Therefore, it 
can be stated with confidence that FAHT’s nonlinear boundary condition routine is 
reliable and accurate. 

r J» ble 3 Vi , Comparison of FAHT and Sinda-87 Temperature Predictions, 

(x,t) - I at I ime = 2 seconds (Radiation Boundary Condition Verification). 


Code 

x = 0 mm 

0.625 

1.25 

1.875 

2.5 

FAHT 

3002 

2410 

1887 

1448 

1098 

SINDA-87 

3001 

2411 

1888 

1449 

1099 


- Test problem with T gas = 3600 K, e = 0.8 

- 320 elements, At = 0.1 sec 
Constant properties 
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III. 3 Variable Propety Routine 


The radiation test case is extended to one with a variable thermal conductivity with 
k = 3.98 + 0.002 T, while density and specific heat are kept constant. The specified 
thermal conductivity function will result in an increase in k and thermal diffusivity 
by a factor 2.44 over the temperature range of 300 to 3600 k. It should be added that 
it is not necessary to have all the thermophysical properties as variables in order to 
verify the variable property routine. This is due to the fact that the same logic and 
routine is used to evaluate the new property values. 

In section III. 1 Eq. (9) was proposed for selection of a time step (for a given Ax) which 
would ensure a realistic answer, and upon further mesh refinement would lead to an 
accurate solution. Thermal diffusivity is a parameter in Eq. (9). For variable 
property problems, it is recommended to use the maximum and minimum values of a 
to obtain the corresponding limits on At. The maximum a value will result in a 
smaller allowable time step range. However, it should be noted that the thermal 
penetration front moves into the undisturbed domain at a temperature which is less 
than the maximum value. Once a region is penetrated with the thermal front, the 
finite element solution routine is not prone to result in unrealistic temperature 
values due to subsequent changes in a. Therefore, the value of a at the highest 
temperature of the region penetrated during a given At is the controlling parameter. 
Usually this controlling temperature is much lower than the maximum value. 

Based on k = 3.98 W/m-K and Ax = 0.625mm Eq. (9) results in 0.07 < At < 0.21 
sec. The highest possible value of thermal conductivity (atT = 3600 k) is 11.18 
W/m-K. The allowable At range corresponding to this k value is 0.025 < At < 

0.075 sec. As discussed earlier the controlling temperature is indeed much smaller 
than the maximum value (3600K). Therefore, At < 0.075 sec, is not an accurate 
estimation of the upper limit on At. Thus, the variable property test problem 
solution is obtained for Ax = 0.625mm and At — 0.1 sec. (same as previous cases), 
and the results are tabulated in Table 4 along with SINDA-87 prediction. 

The maximum variation between FAHT and SINDA-87 results is less than 0.8%. 

I hus, it can be stated with confidence that FAHT’s variable property routine is 
reliable and accurate. 
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Table 4. Comparison of FAHT and SINDA - 87 Temperature predictions, 
T (x,t) - T, at Time = 2 seconds (Variable Property Verification). 


Code 

x — 0 mm 

0.625 

1.25 

1.875 

2.5 

FAHT 

2691 

2341 

1997 

1667 

1362 

SINDA 87 

2705 

2357 

2012 

1371 

1680 


Test problem with T . = 3600 K, e = 0.8 

- Thermal conductivity variable, k = 3.98 + 0.002 T 
constant density and specific heat 

320 elements, At = 0.1 sec 
111.4 Application to MNASA Nozzle 

The variable property and nonlinear boundary condition routines of FAHT have 
been verified. Moreover, a procedure is proposed for obtaining a realistic transient 
solution from FAHT. It is further shown that with mesh refinement the procedure 
will converge to an accurate solution. 


In this section, a course mesh model of MNASA nozzle will be used to obtain the 
temperature solution for variable property case with convection and radiation 
boundary condition. Since the results cannot be compared with any reliable 
solution, no attempts will be made at mesh refinement for obtaining an accurate 
so.ution. 1 he course mesh model of MNASA nozzle along with the description of its 
various materials is depicted in Fig. 5. 

1 he largest element size in the radial direction along the exposed surface to the 
thermal load is 0.547 in. Thermal diffusivity of all the materials (i.e., carbon 
phenolic, glass phenolic, silica phenolic, NBR rubber and steel) maybe relevant in 
the calculation of the allowable time step via Eq. (9). If the transient time is such 
that the thermal front reaches an interior material, then a of that region should be 
considered. Therefore, each material has its own restriction on the allowable time 
step. Moreover, the element size in the direction of heat flow may change from 
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Figure 5. Coarse Mesh Model of MNASA Nozzl 


material to material as well as within a given region. Variation of a with 
temperature is another factor which complicates the selection of proper time step. It 
is not in the scope of this study to calculate all the potential applicable ranges for the 
allowable time step. Therefore, the allowable time ranges for carbon phenolic and 
glass phenolic regions based on a values at about lOOCPR and Ax = 0.547 in are 
calculated and presented. At range based on Eq. (9): 

Carbon Phenolic: 90 < At < 270 sec; Glass phenolic : 770 < At < 2300 sec. 

a = 1-1 x 1 0 in 2 /sec; a = 1.3 x 10 4 in 2 /'sec 

The criterion for carbon phenolic indicates that At > 90 sec. However, an estimate 
of location of the thermal front should be obtained to see if the glass phenolic region 
will be effected for time > 90 sec. The estimated penetraction depth for a = 1.1 x 
10 :i in 2 /sec (carbon phenolic) and At = 120 sec is 1.26 in. It should be emphasized 
thato — 1.26 in is an estimate of the location of the thermal front. This calculation 

shows that the thermal front at time =120 sec will be very close to the glass 
phenolic region. Therefore, the thermal diffusivity of this region may be relevant. 

A FAHT run for the model with At = 120 sec resulted in a number of temperatures 
below the initial condition value. This indicates that the thermal front does indeed 
reach the glass phenolic region and the At restriction for this region should be 
considered (7/0 <[ At < 2300 sec). It should be added that the At restriction of this 
region is an estimate based on element thickness of 0.547 in and a at 1000°R. The 
actual element thickness in this region is less than 0.547 in, and a At < 770 sec will 
be acceptable. As stated earlier, the objective is to show the details of the procedure 
lor At selection rather than tedious calculation. 

FAHT’s predicted temperature field for the coarse mesh model of MNASA nozzle 
with At = 700 sec after one time step is shown in Fig. 6. The predicted temperature 
field seems i easonable. I n a bsence of a reliable solution to compare with no other 
conclusion can be drawn. The maximum predicted temperature at one node is 
6007"F which is greater than T r = 6000"R. This is obviously an error. This error 
may be due to the input mesh geometry. Detail study of the mesh has revealed that 
i n some i egions the nodes of adjacent elements do not coincide. Moreover, there are 
elements which are triangular elements with four distinct nodes. FAHT can not 
account for traingular elements when the third and fourth nodes do not coincide. 

I n fortunately , the limited duration of this project did not allow for mesh correction 
and further investigation of the source of this error. 
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IV. CONCLUSIONS AND RECOMMENDATIONS 

I V . 1 Conclusions 

The heat transfer module of FANTASTIC code (FAHT) is studied and evaluated to 

the extend possible during the ten weeks duration of this project. The conclusions of 

this work are: 

• It is established that with improper choice of element size and time step FAHT’s 
temperature prediction at some nodes, will be below the initial condition value. 
The source of this unrealistic temperature prediction is identified and a 
procedure is proposed for avoiding this phenomenon. It is further shown that the 
proposed procedure will converge to an accurate prediction upon mesh 
refinement. 

• Radiation boundary condition solution routine of FAHT is verified. 

• Variable property solution routine of FAHT is verified. 

• Verification of the ability of FAHT to model convection heat transfer in a porous 
domain, independent of pyrolysis process , is not possible. 

• FAHT users are advised to bypass the transient logic of FAHT by specifying a 
fixed time step based on the proposed criteria. 

• Experienced and dedicated personnel working as a team are required for 
successful usage of FANTASTIC Code. 

• 1 be temperature field solution will have some errors due to lack of explicit 
treatment of the convection term. 

• The pore pressure equation, as given by Eq. (3) is nonlinear. This equation, 
however, is treated as a linear equation in FAHT. Therefore, the pore pressure 
solution will not be accurate, which will result in an inaccurate velocity solution. 

I he inaccuracy associated with the velocity solution, in turn, will increase the 
error in the temperature solution. 

• The How of pyrolysis gases through the char zone are modeled as a non-reacting 
low (frozen). Therefore, the endothermic reactions which take place in this zone 
can not be accounted for. Furthermore, the changes in th porosity and 
permeability in this zone cannot be calculated and accounted for. 

• The permeability of the material is assumed to have the same value in the virgin, 
decomposition and char zones. As mentioned earlier, the permeability is 

changing rapidly with time in the decomposition zone. Nonetheless, FAHT can 
not account for this variation. 
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• There are no provisions for accounting for the initial porosity of the virgin 
material. 

• The momentum equation is based on the Darcy law. However, it is known 
(NASA CR-1903, 1971 ) that the inertial effects play an important role due to 
relatively high mass fluxes of the degradation products. An accurate modeling of 
the momentum equation requires the inclusion of the Forchheimer term. 

Another motivation for the use of Forchheimer-extended Darcy equation of 
motion for flow through porous media is the following. One of the objectives of 
the exploratory test program of the Solid Propulsion Integrity Program (SPIP) is 
to provide empirical relations for permeability of various candidate materials. 
There may be cases where the permeability should be determined via the 
Forchheimer-extended Darcy equation. In that case we do not have any analysis 
tool that can properly use the permeability data. 

IV. 2 Recommendations 

The recommendations of this study have two objectives. First, the issues which 
should be worked in order to get FAHT in an operational status. The second 
objective has a long term view of the required capabilities for accurate analysis of 
charring-decomposing ablator materials including future candidate composites. 

Immediate issues which should be addressed in order to get FAHT in an operational 
status are: 

• Fixing of the identified bugs. 

• Details of the modeling and programming of the convection heat transfer of the 
pyrolysis gas as well as the mass diffusion equation should be provided by FaAA 
for accuracy analysis. 

• Verification of pyrolysis modeling (without ablation) through a comparison of a 
one-dimensional problem with CMA. 

• Verification of FAHT’s ability to model the moving ablating surface through a 
comparison of a one-dimensional problem with CMA. 

• Test run of an MNASA motor model with a refined mesh. 
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I he motivation for recommending a long-term plan of work is as follows. The 
problem of analysis of a charring-decomposing material with an ablating surface is 
very complex and difficult. Experimental verification of any code is difficult, 
expensive and at best will provide a very rough, within the range, comparison. At 
the present time we do not have a " research ” tool which models the energy and pore 
pressure equations accurately without any expedient simplifications for ease of 


programming. Moreover, we do not have a tool which models a dynamic permeability 
and can account for reactions between the pyrolysis gases and the char zone. 
Additional concern is that we do not a code whose momentum equation includes the 
inertia! effects which can play an important role when the mass fluxes are relatively 
high. In short, we do not have a research tool that includes all the known effects. 
Therefore, the research, development, design and testing verification can not 


proceed in a systematic manner. A research code, not a user friendly code or 
necessarily computationally efficient code, is needed to approach the problem in a 
systematic fashion. In that case, we can establish and distinguish primary, 
secondary and tertiary effects. This research code should have provisions to 


accommodate the data about the behavior of new candidate materials as becomes 
available. Concurrent with development of this research code, experimental work 


must proceed to provide extensive data about thermophysical properties, 
permeability and the possible reactions that take place in the char zone. 
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